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ABSTRACT 

Observational evidences point to a common explosion mechanism of Type la supernovae based 
on a delayed detonation of a white dwarf. However, all attempts to find a convincing ignition 
mechanism based on a delayed detonation in a destabilized, expanding, white dwarf have been 
elusive so far. One of the possibilities that has been invoked is that an inefficient deflagration 
leads to pulsation of a Chandrasekhar-mass white dwarf, followed by formation of an accretion 
shock that confines a carbon-oxygen rich core, while transforming the kinetic energy of the 
collapsing halo into thermal energy of the core, until an inward moving detonation is formed. 
This chain of events has been termed Pulsating Reverse Detonation (PRD). In this work we 
present three dimensional numerical simulations of PRD models from the time of detonation 
initiation up to homologous expansion. Different models characterized by the amount of mass 
burned during the deflagration phase, Afdcfl, give explosions spanning a range of kinetic energies, 
K - (1.0 - 1.2) X 10^1 erg, and ^^Ni masses, M{^'^m) - 0.6 - 0.8 M©, which are compatible 
with what is expected for typical Type la supernovae. Spectra and light curves of angle-averaged 
spherically symmetric versions of the PRD models are discussed. Type la supernova spectra pose 
the most stringent requirements on PRD models. 

Subject headings: Supernovae: general - hydrodynamics - nuclear reactions, nucleosynthesis, abundances 



1. Introduction 

Type la supernovae (SNIa) are the most en- 
ergetic transient phenomena in the Universe dis- 
playing most of its energy in the optical band 
of the electromagnetic spectrum, where they ri- 
val in brightness with their host galaxies during 
several weeks. The importance of SNIa in as- 
trophysics and cosmology is highlighted by their 
use as standard (or, better, calibrable) candles to 
measure cos mic distances an d related cosmological 
parameters dPerlmutter e~al. .l998; Schmid t et al 



19981 : iRiess et al.l 120011 ). However, in order to 



achieve a high precision in the distance determi- 
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nation as required, for example, to determine the 
equation of state of the dark energy component 
of our Universe, it is necessary to understand the 
physics of SNIa explosions. From the theoreti- 
cal point of view, the accepted model of SNIa 
consists of a carbon-oxygen white dwarf (WD) 
near the Chandrasekhar mass that accretes mat- 
ter from a companion in a close binary system. 
This model accounts for the SNIa sample homo- 
geneity, the lack of prominent hydrogen lines in 
their spectra, and its detection in elliptical galax- 
ies. Moreover, massive WDs are extremely un- 
stable bodies in which a modest release of en- 
ergy can produce a huge change of size (i.e. the 
ratio of binding to gravitational energy is only 
~ 15%, while in a normal star it is ^ 50%). There 
are two main ingredients of the standard model 
that are still poorly known: the precise configura- 
tion of the stellar binary system and its evolution 



prior t o thermal runaway o f the WD (ILanger et al 



2OOOI : iNomoto et aD 120031 : IPiersanti et al.l l2003t 



Han fc Podsiadlowskil 12001 iBadenes etaU 12007 : 
Hachisu et al."2008'), and the explosion mechanism 
(^Hillebrandt & Niemcyer 2000). Fortunately, as 
long as a carbon-oxygen WD reaches the Chan- 
drasekhar mass, its previous evolution is not crit- 
ical for the explosion because the WD structure 
is determined by the state of a degenerate gas of 
fermions that can be described with a single pa- 
rameter: the mass of the star. This fact leaves 
the explosion mechanism as the most relevant un- 
known concerning SNIa. 

In spite of continued theoretical efforts dedi- 
cated to understand the mechanism behind SNIa, 
realistic simulations are still unable to provide 
a satisfactory description of the details of these 
thermonuclear explosions. Nowadays, there is 
consensus that the initial phases of the explo- 
sion involve a subsonic thermonuclear flame (de- 
flagration) , whose propagation competes with the 
expansion of the WD. After a while, the cor- 
rugation of the flame front induced by hydro- 
dynamic instabilities culminates in an accelera- 
tion of the effective combustion rate. However, 
the huge difference of lengthscales between the 
white dwarf 10^ cm) and the flame width 
1 cm) prevents that large-scale numerical sim- 
ulations can resolve the deflagration front, mak- 
ing it necessary to implement a model of the sub- 
sonic flame. The precise value of the effective 
combustion rate is currently under debate, as it 
depends on details of the flame model. Recent 
three-dimensional (3D) models calculated by dif- 
ferent groups have shown that pure deflagrations 
always give final kinetic energies that fall short 
of 10^^ ergs, while leaving too much unburn t car- 



bon and oxygen close to the center (jR cineck e et al 



2002 : 
2005; 



Gamezo et al.l 



Ropke et al.l 



2003t 
200<fl . 
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Because both signa- 
tures are at o d ds with observati onal constraints, 
Gamezo et aD (|200l l2004l |2005[ ) concluded that 
the only way to reconcile three-dimensional simu- 
lations with observations is to assume that a det- 
onation ignites after a few tenths of a solar mass 
have been incinerated subsonically (delayed deto- 
nation). 

Through a systematic an alysis of a large sam - 
plc of weU observed SNIa, iMazzah et "all (|2007l ) 



showed that in spite of the significant dispersion of 
supernova luminosities, implying different ejected 
masses of ^^Ni, other characteristic features of 



SNIa are remarkably homogeneous. In particular, 
the maximum velocity reached by intermediate 
mass elements, 11000 kms^^, indicates a sim- 
ilar extent of thermonuclear burning in all SNIa. 
Admitting that part of the diversity shown by the 
SNIa sample is directly related to the explosion 
mechanism of Chan drasekhar-mass white dwarfs 
( Hatano et aLllioOoh . multidimensional hydrody- 
namic simulations of the phenomenon should be 
able to explain the wide range of nickel masses 
that are seemingly synthesized in different events. 
Lacking multidimensional calculations of the light 
curve and spectra of these models, it is certainly 
risky to compare directly the outcome of three di- 
mensional explosion simulations with SNIa obser- 
vations, in order to elucidate the more promising 
explosion mechanism and discriminate between 
models. However, even acknowledging that a per- 
fect match between the observed luminosity and 
spectral evolution of SNIa and those deduced from 
angle-averaged versions of the three dimensional 
explosion models is not to be expected, there are a 
few elementary requirements that a model of typ- 
ical SNIa should pass: 

• There should be no more than a few tenths 
of unburnt carbon moving at low veloci- 



ties (iMarion et al.l l2003t iBaron et all 120031: 



Kozma et al. 200 



The upper limit of the ejected mass of ^^Ni 
should reach ~ 1 Mq, after neutronization 
due to electr on captures is prop erly taken 
into account ( Mazzali et al. 2003). 



Intermediate-mass elements should be gen- 
erated in amounts exc eeding ~ 0.2 Mq 
(iGarcia-Senz et al.|[2007h . 



• The ejecta should be chemically stratified, 
as demanded by obser vations of supe r novae 



and remnants (e.g . Badenes et al.l 20061 : 



Gerardv et al.ll2007l : iMazzah et al.ll2008f ). 



• Large clumps of radioactive ^^Ni and other 
Fe-group elements should not be present at 
the photos phere at the time o f maximum 
brightness ([Thomas et al.ll2002r ). 



• In order to reproduce the observed width- 
luminosity relationship the kinetic energy 
and the total mass burned have to be in 



2 



the ranges 
and Mburnt 



K 



( Wooslev et al 



~ (1.0- 
1.0 - 1.2 
200710. 



1.4) 

Mo, 



X 10^^ ergs 
respectively 



Baron et al.l ( 2008h calculated detailed non- 
LTE synthetic spectra of angle-averaged spher- 
ically symmetric versions of Pulsating Reverse 
Detonation (PRD) models. As compared with 
early, maximum, and post-maximum spectra of 
well observed SNIa, the PRD models produced an 
excessively prominent C II feature, and the spec- 
tra appeared too red due to the pres ence of iron- 



peak elements in the outer layers. Baron et al 



(j2008l ) pointed out that fitting of real SNIa spectra 
might improve if some fraction of the iron-peak el- 
ements moving at high velocities were substituted 
by intermediate-mass elements. 

In a companion paper (Bravo & Garcfa-Senz, 
2009; hereafter paper I) we have analyzed the 
conditions for detonation ignition after a pulsat- 
ing phase of a white dwarf. A confined detona- 
tion ignition (GDI) can be induced by an accre- 
tion shock, formed around a carbon-oxygen rich 
core after pulsation following subsonic burning of 
a mass Mdefi < 0.3 M©. The mechanism of GDI 
is robust, and detonation ignition conditions are 
reached non-marginally within the shocked flow 
for a wide range of initial configurations. How- 
ever, if the nuclear energy released before the pul- 
sation is close to the white dwarf binding energy, 
-^dcfl ~ 0.29 Mq, the accretion shock is not effi- 
cient enough to initiate a detonation in a mixture 
of carbon and oxygen. In that paper we have also 
shown that if G-0 matter is contaminated by a 
small amount of He (X(He) - 0.05 - 0.10) the 
conditions for a successful GDI are relaxed sub- 
stantially, making even easier the detonation of 
the core. A small cap made of He might rest atop 
of the white dwarf at the moment of thermal run- 
away as a result of previous accretion from the 
companion star in the binary system. This he- 
lium could cither have been accreted directly from 
a He star or result f rom the nuclear pro c essin, 
of accreted hydrogen. iHoflich fc Khokhlovl (199 
suggested that the presence of a ~ 0.01 cap of 



^Wooslev et al.l | |2007| ) also required that the inner layers of 
ejecta (< 0.8 Mq) should be fairly well mixed, but such 
extent of mixing can be in conflict with what is required to 
explain SNIa spectra both in the optical and in the infrared 
dMotohara et al.ll200d iGerardv et alll2007l : iMazzali et"al] 

l2008h 



He would improve the match of calculated and ob- 
served light curves. There have been some claims 
of detection of He I lines at 2.04 yum and 1.052 ^m 
in the spectra of the spe ctroscopically normal S N 
1994D and o ther SN Ia (iGumming et al.lll994bl la: 
Nomoto et al. 2003; Pignata et al.l2008l). although 



its identification is very unc e rtain (IWheeler et aL. 
19981 : iMazzah fc Lucvlll998t iHatano et al.lll999h . 



In a few test calculations shown in the present pa- 
per we have assumed the presence of such a small 
cap of He. 

In the present work we give details of the det- 
onation phase of PRD models computed in three 
dimensions, together with an analysis of the nu- 
cleosynthesis and the bolometric light curve ob- 
tained from angle-averaged spherically symmetric 
versions of such models. We also discuss the conse- 
quences that departures from spherical symmetry 
of the ejecta can have over the spectra. In Sec- 
tion 2 we explain the methodology used to com- 
pute the three-dimensional supernova models pre- 
sented in this paper. In Section 3 the evolution 
of the white dwarf during the detonation phase is 
described, while in Section 4 the properties of the 
homologous expanding ejecta are analyzed. Fi- 
nally, our conclusions are given in Section 5. 

2. Methods and models 

The simulations were performed with a three- 
dimensional Lagrangian Smoothed Particle Hy- 
drodynamics (SPH) cod e suited to the modelling 



1998 



of Type la supernovae (|Garcia-Senz et al. 
iGarcia-Senz et all Il999l ) . Gontrary to Eulerian 
methods in which usually the spatial resolution 
is a free numerical parameter and the mass reso- 
lution is determined by the local density, in SPH 
it is the mass resolution of the calculation which 
IS fixecH. In SPH, the number of interpolating 
neighbors is what determines the spatial resolu- 
tion of the calculation. In our models, the num- 
ber of neighbors was variable in the range 40-90, 
while the total number of equal-mass particles was 
N = 250, 000. The maximum spatial resolution 
reached during the simulations of the detonation 



^Although it is possible to dynamically increase the mass 
resolution during the SPH calculation using a technique 
known as particle splitting, that is the equivalent to AMR 
methods in Eulerian codes, we have not used it in the cal- 
culations reported in the present paper. 
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phase of the PRD models was 9—12 km, depending 
on the model (the radius of the star at that epoch 
was ~ 40, 000-95, 000 km). Further details of the 
numerical techniques applied in the SPH code are 
provided in the appendix. 

The physical ingredients included in the hy- 
drocode allow to describe accurately all the phe- 
nomena relevant to thermonuclear supernovae. 
Gravity was calculat ed by means of the hierar- 
chical tree method ( Hernquis^ 1987 ). retaining 
up to the quadrupole terms in the multipolar ex- 
pansion. In order to improve the accuracy of 
the calculation, the gravitational interaction be- 
tween close particles was computed directly. The 
equation of state includes contributions of par- 
tially degenerate and relativistic electrons with 
pair corrections, ions considered as an ideal gas 
plus Coulomb corrections and radiation. The nu- 
clear binding ener g y and electron capture rates 



([Fuller et all Il982t iMartfnez-Pinedo et al.l l2000l ) 



on matter in nuclear statistical equilibrium (NSE) 
were interpolated from a table that uses density, 
temperature and electron mole number as input. 
Weak interactions in general, and electron cap- 
tures in particular, are necessary for an accurate 
description of thermonuclear supernovae because 
they affect the dynamical evolution by decreas- 
ing the total number of electrons that support 
the white dwarf structure. The mean electron 
mole number of the ashes, Al^^, determined by 
electron captures during the deflagration phase 
is ~ 0.484 mol g~^. This neutronization 

translates into a decrease of electron pressure of 
as much as l^pl/p ^ (8/3) AF^f ~ 4%, of the same 
order as the increase in pressure due to incinera- 
tion of carbon and oxygen at a density of 10^ g 

Cm-3, Apburn/P 10%. 

The nuclear network consumed a huge part of 
the computational resources chiefly due to the very 
small time-steps required to integrate accurately 
the nuclear kinetic equations, and because at some 
point of the calculation there were tens of thou- 
sands of particles undergoing simultaneously nu- 
clear combustion. A small network of 9 isotopes: 
n^ilLiHeJfC, 1^0, ^ONe, ^^Mg, ^Ssi, and ^^^Ni 
(jWooslevI [l98fil) allows to calculate the approxi- 
mate nuclear energy input and a rough nucleosyn- 
thesis with a moderate load of CPU time. The 
pe rformance of th i s netw ork was extensively tested 
in Timmes et al. ( 2000h . An operator-split ap- 



proach was used to make the calculation feasible. 
As soon as the characteristic nuclear time-step of 
a given particle became smaller than a prescribed 
fraction of the dynamical time-step its chemical 
evolution was followed isochorically, and both the 
nuclear network and the equation of state were 
decoupled from the hydrodynamic evolution until 
the dynamical time-step was recovered. When the 
temperature of a mass element reached 5.5 x 10^ K, 
NSE was assumed. Once attained, NSE was main- 
tained as long as the temperature stayed above 
2 X 10^ K, providing a nuclear energy generation 
rate accurate enough for the hydrodynamical sim- 
ulations. Weak interactions during NSE deter- 
mine the electron mole number, Y^, whose evo- 
lution was calculated at each timestep by solving 
the equation: dYc/dt — T^iXiYi, where stands 
for all kind of weak interactions: Ai > for (3~ dis- 
integrations and e+ captures, and Ai < for /3+ 
disintegrations and e~ captures, while the molar 
fractions, F^, were set by the NSE equations. 

2.1. Testing the code 

We have tested the performance of our SPH 
hydrocode simulating two three-dimensional de- 
layed detonation scenarios available in the scien- 
tific literature, which had been computed using 
quite different numerical methods (both for resolv- 
ing the hydrodynamical equations and to describe 
the nuclear flames). We have chosen delayed det- 
onation models that bear many similarities with 
the PRD scenario, i.e. those in which a delayed 
detonation is formed at a low density and propa- 
gates inwards through a carbon-oxygen rich white 
dwarf core. Our reference models are the Gravita- 
tion ally Confined Detonation (CCD) model YlOO 
from lPlewa ( 2007 ). and the tu rbulent delayed det- 
onati on model DD_005 from ijRdpke fc Niemeyer 
20071 ). The first one begins with the incinera- 
tion of a small off-centered bubble which floats to 
low-density regions consuming a few hundredths 
of a solar mass to finally break-up the white dwarf 
surface. The combustion ashes slide through the 
surface to finally concur at the opposite pole ig- 
niting a detonation which processes most of the 
white dwarf to NSE, mainly in the form of ^^Ni. 
The turbulent delayed detonation model DD_G05 
is based on a physical criteria of detonation igni- 
tion, namely the transition from a laminar flame to 
the distributed regime of nuclear combustion. We 
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have computed models GCDl and TURB7 using 
similar initial models and physical assumptions as 
for YIOCH and DD.005, respectively. Model GCDl 
started from a single incinerated bubble of radius 
53 km located 100 km off-center. Model TURB7 
started from the same initial model as model DF29 
in paper I, but it incorporated a prescription for 
turbulent acceleration of the burning rate accord- 
ing to the local velocity fluctuations. Although the 
initial configuration of TURB7 is not the same as 
in DD_005, both models burnt similar masses sub- 
sonically, that were distributed in a small number 
of ash clumps, and used the same criteria for the 
deflagration-to-detonation transition. As we want 
to test the ability of the code to deal with det- 
onation propagation, the similarity of the global 
properties of the white dwarf at the end of the de- 
flagration era is what makes meaningful the com- 
parison of models DD_005 and TURB7. 

The outcome of the additional models is sum- 
marized in Table [TJ where we give the amount of 
mass burned during the initial deflagration phase 
(when available), the final kinetic energy, and the 
ejected mass of ^^Ni, NSE elements, carbon plus 
oxygen, and intermediate-mass elements (details 
of these simulations and the numerical setup used 
will be presented in a forthcoming paper, Garci'a- 
Senz et. al. 2009). The mass of ^^Ni in the ref- 
erence models is just an upper limit because elec- 
tron captures were not included in the hydrody- 
namic calculations. As can be seen in the Table, 
the matching between our models and the refer- 
ence ones is excellent, especially when one com- 
pares GCDl with YlOO, which are nearly twin 
models. Model TURB7 reproduces quite well the 
main properties of model DD_005, but some small 
differences remain in the synthesized masses of 
intermediate-mass elements (IME) and NSE ele- 
ments. The reason for these discrepancies may 
reside in different approaches to the description 
of the flame as well as on the effect of the elec- 
tron captures on the dynamics of the explosion. 
Even though neither DD_005 nor YlOO included 
electron captures, the amount of mass burned at 
high densities during the deflagration phase was 
much larger in the turbulent delayed detonation 
model (0.38 Mq vs ~ 0.06 M©), hence the impact 



of the neutronization on the explosion dynamics 
should be also larger in TURB7. The final mean 
electron mole number in GCDl and TURB7 was 
= 0.4964 and = 0.4932, respectively, i.e. the 
neutronization in the turbulent model was about 
double than in the other one. 

2.2. Pulsating reverse detonation simula- 
tions 

Our working hypotheses are that during the 
deflagration phase it is burnt an insufficient 
amount of mass to unbind the white dwarf, 
and that the outcome of the explosion is pre- 
dominantly determined by Mdefl. The simula- 
tions of the deflagration phase started from ini- 
tial models consisting on spherically symmetric 
cold (T = 10^ K) isothermal white dwarfs built 
through a fourth order Runge-Kutta integration. 
The one-dimensional profile was mapped to a 
three-dimensional distribution of N particles and 
af terwards relaxed with th e procedure explained 
(|l998h 



Garcia-Senz et al 



Once a stable enough 
configuration was obtained, the temperature of a 
small subset of particles was artificially increased 
up to the self-consistent value of NSE tempera- 
ture at such density (assuming isochoric adiabatic 
burning) and the ensuing dynamical evolution was 
followed by means of the SPH hydrocode. 

The main properties of our PRD models are 
given in Table [H For clarity, we have changed the 
name of the models with respect to previous pub- 
lications. In the present work, PRD models are 
designated with the acronym PRDnn, where 'nn' 
gives the hundredths of a solar mass burned dur- 
ing the initial deflagration phase. For instance, in 
model PRD18 the mass burned subsonically before 
the pulsation was M ^p^ = 0.18 M© (this model 
was called PRD6 in iBravo fc Garcfa-Senj l2006l 



and Baron et al. 20081 while the present mod el 



PRD14 was called PRD5.5 in lBaron et ani2008[ ). 



^Note that YlOO is a two dimensional model, while GCDl 
was computed in 3D 



The evolution of the central density and the 
chemical composition of model PRD14 is shown 
in Fig. [T] At the beginning combustion propa- 
gated subsonically as a deflagration wave, with 
an effective combustion rate low enough to en- 
sure that the flame quenched before the white 
dwarf became unbind. During the first seconds 
of the ensuing pulsation the incinerated bubbles 
fioated to the surface, leading to composition in- 
version, i.e. the internal volume, which we call 
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hereafter the core, was plenty of cold fuel, i.e. car- 
bon and oxygen, while the ashes of the initial com- 
bustion, mostly hot iron and nickel, were scattered 
around (see Fig. [2]). Shortly after, the expansion 
of the core came to an end while part of the enve- 
lope re-collapsed. Eventually, an accretion shock 
was born because of the impact of the in-falling 
material onto the carbon-oxygen core. Upon re- 
collapse nuclear reactions light again a few tenths 
of a second before the time of maximum compres- 
sion (~ 4.4 s). In this paper, our interest is fo- 
cused on the events that ensue after a detona- 
tion is ignited due to the energy transferred to 
the core through the accretion shock (the defla- 
gration phase has been th oroughly described in 
Garcfa-Senz fc Bravdl2005l ). 



3. Description of the detonation phase of 
PRD models 

In paper I we have shown that after a pulsa- 
tion the white dwarf reaches conditions favorable 
for a GDI if Mdcfl is well below 0.3 M© and 
nuclear reactions are turned off. However, when 
nuclear reactions are included in the calculation 
the release of nuclear energy contributes not only 
to temperature increase but also to raising the gas 
pressure, thus disturbing the dynamical evolution 
during the collapsing phase of the pulsation. In 
the present simulations we have allowed the det- 
onation to ignite spontaneously at the time and 
place determined by the SPH code as a result 
of the hydrodynamical evolution coupled to the 
nuclear network. This approach must be looked 
with caution as we cannot resolve the detonation 
front thickness, hence our simulations might not 
describe accurately the process of detonation for- 
mation. However, as we will see later in Sect. 21 
the explosion properties are quite independent of 
the details of the detonation initiation, hence our 
approach is justified. 

3.1. Initial steps of detonation propaga- 
tion 

It is interesting to analyze the initial steps of 
detonation propagation, keeping in mind the cau- 
tionary remark given above about the capabili- 
ties of our code to describe accurately the birth 
of a detonation front. Figure [3] shows the evolu- 
tion of a zoomed slice across model PRD 14 for 



two times: just after the nuclear timescale be- 
comes shorter than the hydrodynamical timescale, 
and 0.1 s later (top images). The maximum tem- 
peratures achieved in the region are high enough 
(T > 5 X 10^ K) for explosive silicon burning, 
whereas carbon and oxygen are destroyed in a 
shorter timescale once the temperature rises above 
~ 3 X 10^ K as highlighted by the fuel contours. 
The thermal gradient in the radial direction is 
large, on the order of 100—200 K cm~^, although it 
is influenced by the resolution of the numerical cal- 
culation. Shock fronts propagate ahead and close 
to the thermal waves, as revealed by the sculpt- 
ing of the isobars following the contour of the high 
temperature regions and by the large pressure gra- 
dient in front of the thermal front (bottom left 
image): the detonation is born. 

The performance of the detonation capture al- 
gorithm is illustrated in the bottom right image 
on Fig. [21 where it is shown the region where the 
adaptive elliptical interpolator gives a sharper det- 
onation front. As can be seen, a detonation is 
detected all around the thermal wave except the 
outermost part characterized by a shallow thermal 
gradient. Figure |H displays a complementary view 
showing the location (in fact, the projection onto 
the yz plane) of the actual SPH particles within 
the detonation front, i.e. those for which the det- 
onation capture condition has been verified. The 
color map in the same Fig. shows the smoothing 
length, whose value at the detonation front is of 
order 30 km. 

3.2. Burning of the white dwarf core 

In the PRD scenario, the last phase of the ex- 
plosion starts when the converging reverse detona- 
tion wave is launched. Once initiated, the detona- 
tion is self-sustained by the nuclear energy release 
from burning of the carbon-oxygen fuel that per- 
vades the hydrostatic core. 

The evolution of the radii of spherical shells of 
given Lagrangian mass is shown in Fig. [5l together 
with the evolution of the detonation front. Be- 
cause the model does not have spherical symme- 
try, the location of the detonation front is only 
approximate. In that figure, the inner and outer 
boundaries of the detonated region have been de- 
fined as the radius of the shells for which at least 
half of the carbon and oxygen remaining after the 
deflagration phase had been processed by the det- 
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onation wave (see Fig. [6] and the next paragraph 
for a complementary view of the evolution of the 
detonation wave) . It can be seen that the inwards 
propagation of the thermonuclear burning is very 
fast, as corresponds to a detonation front. This 
detonation front processes most of the inner core 
of the white dwarf, but it stalls before arriving 
to the very central regions and fails to burn the 
innermost few hundredths of a solar mass. 

Fig. [S] illustrates the evolution of the white 
dwarf core during the propagation of the reverse 
detonation for model PRD14, with a view com- 
plementary to the one in Fig. [5l At the time of 
the first snapshot, t — 3.8 s, the accretion shock is 
about to form. The C+O mass fraction contours 
maintain a high degree of spherical symmetry in 
the core, increasing their value steeply towards 
the center (lighter contours belong to larger C+O 
mass fractions). In the first snapshot, a hot spot 
is visible in the yz plane dX r ^ 2 000 km, which 
is a result of a slight asymmetry in the collapsing 
flow (the temperature and chemical composition 
maps a few tenths of a second before do not show 
any sign indicating the imminent formation of a 
hot spot). The combination of high temperatures 
and abundant fuel culminated in the formation of 
a detonation already visible in the second snap- 
shot. At the same time, the accretion shock had 
heated the fuel in a hot spherical shell located at 
r - 3000 km up to T > 1.5 x 10^ K. The lateral 
spreading of burning was not due to detonation 
propagation, but it was caused by the accretion 
shock, that continuously fed mechanical energy 
into the core sparking detonations all around the 
hot shell. During the first 0.3 s after detonation 
ignition the accretion shock remained stationary, 
close to the core. Afterwards, the overpressure 
generated by the nuclear energy released pushed 
out the accretion shock, that detached from the 
dense core (last snapshot). In the outer layers of 
the core the density was low enough to allow for in- 
complete silicon burning and leave a composition 
rich in intermediate mass elements, mainly Si, S, 
Ca and Ar, while in the inner regions the burning 
proceeded all the way up to ^^Ni (the central den- 
sity at the moment of formation of the detonation 
was in the range (0.5 — 4) x 10* g cm~'^). 

In our simulations the detonation propagated 
inwards burning all the fuel except for a small vol- 
ume close to the center of the white dwarf. The 



unburnt mass, ^ 0.08 M©, and the central density, 
Pc ~ 4x10® gcm~'^, give a volume of ~ 4x10* km'^ 
or a typical lengthscale of the unburned spot of 
~ 450 km. Our resolution at the center at that 
time was of order 11 km, implying that the det- 
onation failure at the center of the white dwarf 
was not due to insufficient numerical resolution. 
The physical reason for detonation failure at the 
center is the expansion of the external, already 
detonated, core matter. As the detonation wave 
is accessible to the shocked matter through sound 
waves, the rarefaction wave originated by expan- 
sion of the core caught and weakened the shock 
and finally quenched the burning close to the cen- 
ter. It took a few tenths of a second since detona- 
tion initiation to release enough nuclear energy to 
overcome the impact pressure of accreting matter 
and allow the detonated matter to expand appre- 
ciably. As a result of such delay most of the core 
matter was detonated, leaving only a small mass of 
unburned fuel near the center. The amount of un- 
burned carbon and oxygen is sensitive to the den- 
sity structure at the moment of detonation forma- 
tion, being the largest for the less energetic model 
PRD26 (see Table E]). 

3.3. Sensitivity and resolution 

We have checked the sensitivity of the results of 
our simulations to the numerical method, in par- 
ticular to the implementation of artificial viscosity, 
and to the spatial resolution. As it is well known, 
artificial viscosity introduces numerical dissipation 
in the evolution equations, which can generate 
enough entropy at shocks to induce an artificial 
detonation. We have computed two additional 
models starting from Mdefl = 0.18 M©, one with 
the artificial viscosity parameters enlarged by a 
factor two, the other with the parameters divided 
by the same factor. In both cases, the code in- 
cluded an algorithm to inhibit combustion inside 
the artificially widened shock fronts. The detona- 
tion developed in both cases in the same way as 
in our standard calculation, irrespectively of the 
values taken by the artificial viscosity parameters. 
The final kinetic energies and ^^Ni masses synthe- 
sized differred by less than 8% with respect to our 
standard model PRD18. 

We also checked the sensitivity of the results to 
the numerical resolution by running model PRD18 
with the same number of particles but reducing 



7 



the number of neighbors, uq. The smaller number 
of neighbors translated into a smaller smoothing 
length, and thus in an increased capability of the 
code to represent small scale features. Such proce- 
dure allowed us to reach a maximum resolution of 
4 km during detonation propagation. Once more, 
the detonation developed in the same way as in our 
standard calculation. As a drawback, the small 
number of neighbors originated numerical insta- 
bilities in the outermost low-density regions of the 
ejecta. These instabilities prevented the model to 
reach the homologous expansion phase, making it 
impossible to compare the final outcome with that 
of our standard PRD18 model. 

Finally, we checked if the outcome of the ex- 
plosion was sensitive to the presence of a small 
mass of He atop of the white dwarf. This helium 
mixed with the underlying layers of carbon and 
oxygen during the pulsation and was present at 
the detonated layer at ignition time. By adding 
either a total He mass of 0.005 Mq (i.e. a 0.36%) 
or 0.010 Mq (i.e. a 0.72%) to model PRD18 the 
kinetic energy and the mass of ^^Ni synthesized 
in the explosion changed negligibly (< 4%) with 
respect to the standard calculation of the same 
model. Thus, we conclude that the eventual pres- 
ence of He does not introduce additional uncer- 
tainties in the outcome of the PRD models. 

4. Nucleosynthesis, light curves and spec- 
tra: asymmetries, inhomogeneities 

In Table [2] there are given the properties of 
the PRD models we have computed, covering 
the range of subsonically burned masses Mdefl = 
0.14 - 0.26 Mq. The kinetic energy at infinity 
spans the range (1.0—1.2) x 10^^ erg, in good 
agreement with requirements deduced from the 
SNIa light curve shape (see Sect.[T]). In principle, 
it should be possible to obtain PRD models with 
still smaller values of Mdcfi, down to ^ 0.10 Mq or 
even less, which according to the trend displayed 
in Table m would translate into larger kinetic en- 
ergies and ^®Ni masses. However, as the accre- 
tion shock in model PRD14 already forms quite 
close to the total mass of the white dwarf, it is 
to be expected that the outcome of the explosion 
for smaller Afdcfl will not be much different from 
that of model PRD14. On the other side, values of 
Mdefl sHghtly larger than those in model PRD26 



would probably give rise to explosions with slightly 
smaller kinetic energies, substantially smaller '^^Ni 
masses and much larger intermediate-mass ele- 
ments yields. Unfortunately, due to the failure 
to maintain a steady detonation in model PRD29 
we have not been able to verify this last point. 

As seen in the last snapshot of Fig. [HI the det- 
onation smooths the thermal and chemical struc- 
ture of the white dwarf, wiping out the imprints 
left by the deflagration phase in the central vol- 
ume of the ejecta. The detonation also endows 
the mechanical structure of the ejecta with a large 
degree of spherical symmetry. Fig. [7] shows the 
homologous density profile of model PRD14 at a 
time of 34 s since the initial thermal runaway. The 
solid line represents the angle-averaged density ob- 
tained in 1 000 shells evenly distributed in mass, 
while the error bars give the standard deviation of 
the particles density. The small dispersion is an 
indication of the high degree of spherical symme- 
try. The high-velocity low-density external mate- 
rial is best approximated by an exponential den- 
sity law, whereas the central region can be ap- 
proximately described as a nearly uniform density 
core whose mass is ~ 0.82 Mq. The exponential 
tail begins at a velocity ^ 10 000 km s~^, where 
the intermediate-mass elements (Si, S, Ca, ...) are 
abundant: Xime ^ 0.2. The innermost 0.08 Mq, 
composed mainly by unburned C-0, is denser than 
the uniform core by up to an order of magnitude. 

A comparison of the density profiles of the last 
computed models (normalized by t^) is presented 
in Fig. [51 The larger kinetic energy of model 
PRD14 with respect to PRD18 produces a more 
extended tail in the former. The beginning of the 
exponential tail matches quite approximately the 
location where the detonation was ignited in each 
model (see values of Udoto in Table[5]). Whereas the 
transition from the central uniform density region 
to the exponential tail is quite smooth in model 
PRD14, it becomes increasingly abrupt as Mdcto 
and Vdeto decrease. In model PRD26 there is not 
a uniform density core, which is substituted by a 
density inversion below Mdcto- In models PRD18 
and PRD26 the density of the innermost region 
rich in fuel is ^ 4 times larger than the density of 
the nearly homogeneous core, a factor much lower 
than that found in model PRD14. 
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4.1. Nucleosynthesis 

Figs. [9] and \W\ depict three-dimensional views 
of the chemical structure of the ejecta. In Fig.[n]a 
rendering of the final distribution of ^^Ni in model 
PRD14 is shown. The central volume is homoge- 
neously filled with ^^Ni, except for a small ball 
made of carbon and oxygen that is not visible in 
this image, surrounded by intermediate-mass ele- 
ments and unburned carbon and oxygen. Small 
and medium size clumps of ^^Ni and other Fe- 
group elements pervade the outer layers of the 
ejecta. These clumps have a double origin: some of 
them were produced during the initial deflagration 
phase and migrated to the external layers of the 
white dwarf during the expanding phase of the pul- 
sation, while others detached from the incinerated 
core after the end of the detonation phase. These 
clumps are the end product of the impact of the 
detonated core into the surrounding low-velocity 
expanding layers. Thus the final distribution of 
^^Ni comes from two different regimes: a smooth 
inner region, rich in ^^Ni and other Fe-group ele- 
ments, and a clumpy outer region where radioac- 
tive and stable Fe-group and intermediate-mass el- 
ements share the volume with carbon and oxygen. 
A complementary picture of the chemical struc- 
ture is given in Fig. IIOI in which the mean molar 
weight in an octant of model PRD14 is displayed. 
It shows the same gross features seen in Fig. [9l 
but here more details are visible, as for instance 
the central ball made of carbon and oxygen. Nev- 
ertheless, carbon and oxygen are predominantly 
located at the outermost layers of the ejecta. 

In order to compute even coarse light curves 
and spectra of SNIa models and to perform a pre- 
liminary comparison with observations, it is nec- 
essary to obtain the isotopic yields up to yl ~ 60. 
As our SPH code implements just a 9-isotopes nu- 
clear network, we have post-processed the thermo- 
dynamic history of the particles to derive the nu- 
cleosynthetic output of each explosion model. For 
each computed model, the SPH particles were clas- 
sified into three groups: 1) those whose tempera- 
tures never exceeded 10^ K, that were assumed to 
keep their initial composition unchanged, 2) those 
that reached temperatures in excess of Tnse = 
5.5 X 10^ K, that were assumed to attain NSE 
composition and, 3) the remaining particles, that 
were candidates to experience incomplete explo- 
sive burning. The nucleosynthesis of each particle 



belonging to the last group was computed by inte- 
gration of the nuclear kinetic equations along the 
thermodynamic trajectory recorded for that par- 
ticle during the SPH simulations. The method 
followed for the integr ation of the nuclear equa- 
tions was described in ICabezon et al.l (|2004r ). Fi- 
nally, the particles belonging to the NSE group 
were again divided in two categories according to 
the density at which they left NSE. Those par- 
ticles that cooled below 2 x 10® K at densities 
in excess of 10* g cm~^were assigned the chem- 
ical composition in NSE at these temperature and 
density, and at the electron mole number given 
by the SPH model. Conversely, particles that left 
NSE at densities below 10* g cm~'^were assumed 
to experience an alpha-rich freeze-out, and their 
final chemical composition was obtained by inte- 
grating the nuclear network following the cooling 
path of the particle, starting from T = Tnse and 
the chemical composition of NSE matter at that 
temperature. 

The mass of ^^Ni synthesized in each PRD 
model can be found in Table O while Table [3] 
gives the final abundances after radioactive de- 
cays, where we have summed up the abundances of 
the isotopes of each element and we show only the 
most abundant species in the ejecta. The amount 
of ^^Ni produced is in the range 0.6 — 0.8 Mq, 
which is enough to explain the light curves of typ- 
ical SNIa. As pointed before, models with A/dcfl in 
excess of ^ 0.27 Mq might produce lower amounts 
of ^^Ni that might account for sub-luminous SNIa, 
but we have not been able to verify it in the present 
work. The total mass of IME synthesized in the 
three computed models spans a similarly narrow 
range, M(IME) - 0.20 - 0.26 Mq. 

Figures [TTII131 show the final (after radioactive 
decays) chemical profiles of the three PRD models 
as a function of the ejecta velocity. The location 
of the accretion shock, Wacc, and the GDI, Wdcto, 
left no imprint on the chemical profile. The dis- 
persion of the mass fractions agrees with the qual- 
itative view obtained with Figs. [9l and [TOl the in- 
ner volume is chemically stratified while the outer 
layers display large scatters. The transition be- 
tween both regimes takes place at an ejecta ve- 
locity of w ~ 10 000 km s-i in model PRD14 
and V ~ 8 000 km s'^ in model PRD18, while 
in model PRD26 the scatter is already high at a 
few 1 000 km s~^. The abundance of intermediate- 
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mass elements increases towards the surface of the 
ejecta in all models. In model PRD14 the IME 
amount to a significant mass fraction at velocities 
in excess of w ~ 10 000 km s^^, and are dominant 
above v ~ 20 000 km s~^. The large mass fraction 
of Fe-group elements at high velocities is a concern 
because its expected spectral signatures are not 
seen in SNIa. The same is true for models PRD18 
and PRD26. In model PRD18 the abundance of 
carbon plus oxygen increases monotonically up to 
the surface where it amounts to X{G + O) ~ 0.30. 
The C+0 profile in model PRD26 shows a peak at 
the location of the GDI, Vdeto = 8 000 km s~^, and 
declines at larger velocities. Although this model 
has the largest mass of IME, Fe-group elements 
are still dominant through the outermost layers. 

In Fig.[T3]it is shown the detailed chemical com- 
position of model PRD14 after radioactive decays 
as a function of the velocity. We stress that the 
poor chemical stratification of the ejecta is mainly 
due to the mixing between the detonated core ma- 
terial and the surrounding matter that is produced 
during the initial phases of expansion following 
the end of the detonation era. Such a mixing is 
expected to happen in any model in which a det- 
onation propagates inwards from a point close to 
the surface of the star, as is the case in the two 
test models presented in Sect. 12.11 This is illus- 
trated in Fig. [m that shows the chemical profiles 
of model TURB7 at two times: just after the end 
of detonation and when the ejecta is in the ho- 
mologous expansion phase. It is plain to see that 
the chemical species are redistributed through the 
outer ^ 0.6 M© between both times, with the re- 
sult that Fe-group elements replace C-0 and IME 
in the outer layers. Thus, the final chemical pro- 
file of TURB7 displays also a large mass of Fe-Ni 
moving at high velocities. 

4.2. Light curves 

Multidimensional thermonuclear supernova 
models assuming different explosion mechanisms 
have been shown to reproduce reasonably well the 
gross features of broadban d light curves of nor- 
mal Type la supernovae (iBlinnikov et aL 120061: 
Kasen fc Plewal l2007t iRopke et al.l l2007t ). The 
model light curves in these works were obtained 
from one and two-dimensional averaged versions 
of two and three-dimensional supernova models. 
An alternative approach is to compute bolomet- 



ric light curves from spherically symmetric angle- 
averaged versions of three-dimensional SNIa mod- 
els. Bolometric light curves should be less sensi- 
tive to the averaging procedure than broadband 
light curves, as the first ones are not affected 
by shifting of photon wavelengths in and out of 
monochromatic bands, as due for instance to iron 
line blocking. 

Bolometric light curves have been computed 
from angle-averaged versions of the PRD models 
b y means of th e one-dimcnsional code described 
in lBravo et al. I {1993,, 1996). In general, the bolo- 
metric light curves obtained with this code dur- 
ing the pre-maximum and maximum phases are 
in fairly good agreement with those computed by 
dire ctly solving the ra diative transfer equations 
(;see lHoflich et al.lll993[ ). Although our PRD mod- 
els reach the peak of luminosity slightly earlier 
than expected 12 days), the maximum lumi- 
nosity, £max, and ^^Ni mass, M (^^Ni), match al- 
mo st perfectly the empirical relatio nship derived 
bv IStritzinger fc Leibundgu-Q ( 20051 . their Eq. 6). 
In Fig. [16] we show a plot of M(^^Ni) vs the bolo- 
metric light curve decline parameter Amis, both 
from the present three-dimensional models as well 
as from observationa l estimates for a set of well 
observed SNIa (Stritzinger et al.ll2006l ). Models 
PRD14 and PRD26 reproduce nicely the obser- 
vational trend, while model PRD18 declines too 
fast in comparison to SNIa having a similar mass 
of s^Ni. Models GCDl and TURB7 synthesize a 
large amount of ^^Ni and, due to their large kinetic 
energy and the presence of radioactive nuclei up to 
the surface of the ejecta, decline much faster than 
observations seem to indicate (the main parame- 
ter determining the width of the bolometric light 
curve is the photon diffusion time, t^, that is in 
turn aff ected by the kinetic energy: oc K~^/'^, 
rWooslcy"eraI| l2007l ). 



4.3. Spectra 



Baron et al.l (|2008l ) obtained detailed spectra 
of angle-averaged one-dimensional versions of sev- 
eral PRD models, from slightly before up to a 
few days after maximum light. They concluded 
that the synthetic spectra of the PRD models 
did not match the spectra of a few typical well- 
observed SN I a. Th e main disagreements found by 
Baron et al. ( 2008f l were due to the presence of 
Fe-group elements in the outer layers of the ejecta. 
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which translated in an excessive flux at long wave- 
lengths in the computed post-maximum spectra. 
There was also an excess of carbon, that resulted 
in a prominent C II 6580A line at 13 days that was 
too red by ~ 250A and lacked the correct shapcfl. 
In the next paragraphs, we reflect on the impli- 
cations that the asymmetries and inhomogeneities 
of the outer layers of the ejecta might have for the 
predicted spectra. 

Figure [T7] shows a map of the distribution of 
Fe-group elements as would be seen from an ar- 
bitrary line of sight. The column density of Fe- 
group elements in this direction, starting from the 
photosphere, is used to color the map. In order 
to compute the column density, the photosphere 
has been simply deflned as the surface at which 
the optical depth in the line of sight is 2/3, as- 
suming a constant opacity of 0.2 cm^ g~^ (the 
maximum possible column density of Fe-group ele- 
ments down to the photosphere is thus (2/3) /0.2 = 
3.3 g cm~^). This image shows that the distribu- 
tion of Fe through the outermost layers is quite 
asymmetrical, hence their effect on the photon flux 
is expected to be latitude dependent. On average, 
one expects that the asymmetrical distribution of 
Fe should reduce the blocking of UV photons and 
the shifting of the flux to redder wavelengths, es- 
pecially in locations where Fe is scarce such as in 
the upper hemisphere of this image, which should 
allow a better agreement between synthetic PRD 
spectra and observations. 

The implications of the inhomogeneous chemi- 
cal composition for the formation of carbon lines 
are illustrated in Fig. [181 In this flgure we have 
plotted a proxy of the shape of carbon lines from 
different views: the distribution of the absorbing 
mass of carbon between the photosphere and the 
surface as a function of the velocity projected in 
the line of sight (in real supernovae, the depth 
of the absorption lines should correlate with the 
mass of carbon the photons have to cross before 
escaping the ejecta). Obviously, the physics of 
line formation in supernovae is much more com- 
plicated than our simple prescription accounts for, 
but this is enough for the present qualitative dis- 



''We note that the C II line wavelength disagreement might 
be reduced by ~ lOOA if the epoch of the com puted spectra 
were a few days earlier than assumed by Ba ron et a l.!l2008|, 
because the line would form on an outer and faster layer of 
the ejecta 



cussion. What Fig. [T8l shows is that the shape and 
peak of the carbon distribution are quite sensitive 
to the line of sight. The same model is able to 
produce broad and shallow lines and triangular- 
shaped deep lines depending on the line of sight. 
At the same time, the peak of the line can change 
at most by Aw 7000 km s~\ i.e. AA - 150A. 
In spite of the simplicity of these calculations, 
they show that the scattering in the carbon line 
strength due to dimensionality effects might be 
large and, therefore, conclusions obtained by com- 
paring SNIa observations with synthetic spectra of 
angle-averaged versions of three-dimensional mod- 
els have to be taken with care. 

5. Conclusions 

It is amazing that one-parameter spherically 
symmetric explosion models continue provid- 
ing the best fits to SNIa observations. And 
it is even more amazing that elaborate three- 
dimensional simulations of SNIa neither reproduce 
one-dimensional results nor improve on them. In 
spite of recent advances in 3D modeling of delayed 
detonations of white dwarfs, current models are 
still far from passing a suite of basic observational 
tests. The few one-dimensional spectral calcula- 
ti ons from a ngle-averaged 3D SNIa models (e.g. 
iBaron et al.l i2008) did isolate several deficiencies 
when compared to observed SNIa spectra: the 
formation of undetected lines of C, an incorrect 
shape of the S lines, and an excessive abundance 
of Fe-group elements in the outer layers of the 
ejecta. On the other hand, observed bolomet- 
ric light curves are e asily reproduced by current 
3D SNIa models (e.g. iRopke et all 120071 and the 
present work) as they depend mainly on the mass 
of radioactive ^^Ni synthesized during the explo- 
sions. 

In this paper, we have discussed the thermal 
and dynamical evolution of PRD models of SNIa. 
The detonation phase of these models begins af- 
ter accretion shock formation around a carbon- 
oxygen rich core, for a wide range of initial con- 
ditions. The formation of the detonation was ro- 
bustly achieved in our simulations. Even though 
resolution limitations of our large-scale calcula- 
tions preclude us to resolve the detonation wave 
structure, the paths followed by shocked particles 
in the p — T plane surpass clearly the threshold 
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for C-0 detonation, as shown in paper I. Thus, 
a confinement detonation ignition can be safely 
assumed in our large-scale simulations, the excep- 
tion being those models in which the deflagrating 
phase of the explosion, previous to white dwarf 
pulsation, releases an amount of nuclear energy 
close to the binding energy of the star. 

The present PRD models score quite well on 
the basic observational tests quoted in Sect. [T] 

• All models produced less than < 0.13 M0 of 
carbon, and only a small fraction of it was 
ejected with small velocities. 

• The kinetic energy and total burned mass 
stay in the range 1.0 — 1.2 x 10^^ erg and 



1.1 — 1.2 M(n, respective ly (jMazzali et al 



20071 : IWooslev et al.ll2007h 



• The amount of ^^Ni synthesized in the ex- 
plosions ranges from 0.6 Mq to 0.8 Mq that 
is enough to explain normal to bright events, 
although this range appears a little narrow 
to explain the differences in luminosity of the 
whole observational SNIa sample. 

• The mass of intermediate-mass elements is 
larger than 0.2 M0 in all our models. 

• The explosions did not produce big clumps 
of ^®Ni. Instead, radioactive as well as stable 
Fe-group nuclei are assembled in small bub- 
bles pervading the outer layers of the ejecta. 

The chemical profile of the PRD models lacks the 
degree of stratification suggested by observations. 
Although, in general. Si and S display an increas- 
ing mass fraction as function of the ejecta velocity, 
Fe-group elements represent the most abundant 
group except at the innermost layers made of un- 
burnt carbon and oxygen. Stable Fe-group nuclei, 
synthesized at high densities during the initial de- 
flagration phase of the explosion, are not present 
at the center of the ejecta as they should accord- 
ing to infrared observations of SNIa. This might 
be a weak point of the PRD scenario, as bubbles 
floatation and chemical composition inversion at 
the end of the deflagration phase are two of its 
basic ingredients. 

We have pointed out several possible effects 
that the existence of asymmetries in the ejecta 
may have in the supernova spectra. These effects 



might m itigate partia l ly som e of the deficiencies 
found bv lBaron et al. ( 2008f l. who calculated syn- 
thetic spectra using one-dimensional averaged ver- 
sions of the present PRD models. Obviously, much 
work remains to be done before one can ascertain 
if the PRD is a reliable scenario able to account 
for a significant number of SNIa events. 

This work has been partially supported by 
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by the Generalitat de Catalunya 
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A. Numerical techniques for SPH simulations of detonation waves 



Shocked flows present SPH simulations with specific challenges due to the presence of large pressure gra- 
dients and strongly anisotropic distributions of particles. Although there exists extensive literature devoted 
to such topics, we summarize in this appendix the numerical techniques we have implemented in our calcu- 
lations. The interpolating kernel used t o compute the terms in t he 3D SPH evolutionary equations was the 
cubic spline kernel Wij {u, h) defined in lMonaghan fc Lattanzid (|l991l ). 



A.l. Determination of the smoothing length 

We allowed for a variable number of neighbors in order to improve the behavior of the code near shock 
fronts, around which many particles can be crowded leading to fluctuations of the smoothing length if a 
fixed number of neighbors is prescribed. In our implementation, the smoothing length tracked not only the 
number of neighbors but also their distribution around a given particle. The smoothing length, /i^, of particle 
i was defined by: 



no 



N 



(Al) 



where no is a fixed parameter, and rij is the distance between particles i and j. This way, neighbor particles 
contribute to no with a fractional number that is a decreasing function of the distance to particle i, while 
particles farther than 2hi from particle i do not contribute at all to the summation. When neighbor particles 
move across a sphere of radius 2hi the summation is smoothly modified, which translates into a smooth 
adjustment of hi in order for the summation to match no again. Hence, the change of hi is also a sensitive 
function of At and fiuctuations of the smoothing length are avoided. 

In Garcia-Senz et al. we found that a minimum of ^ 40 neighbors in 3D calculations was needed 

to fit the structure of hydrostatic massive white dwarfs without introducing appreciable fluctuations in the 
hydrostatic force. We have found that a value of order no ~ 22 results in a safe minimum of neighbors (> 40) 
that ensures an accurate determination of gradients. Working with a smaller value of no could increase the 
numerical noise beyond acceptable levels, while a signiflcantly larger number of neighbors would smooth 
too much the interpolated variables preventing an accurate representation of the detonation wave. In the 
simulations reported here we have worked with no = 21.875. 

Neighbor searching was implemented through an octal tree, the same used for the calculation of the gravity 
force. Unlike the above described algorithm used to compute the individual smoothing length of each particle, 
which applied the gather method, the evaluation of the kernel in the SPH evolutionary equations followed a 
sy mmetric momentum-cons erving algorithm, combining the gather and scatter methods in the way suggested 
bv lHernquist fc Kat j ( 19891 . see, e.g., their Eq. 2.15). During the PRD simulations reported here the number 
of neighbors contributing to the SPH equations of any given particle lied in the range 40-90. 



A. 2. Artificial viscosity 



Shock waves were handled by artificial viscosity, with the prescription of lBalsaral (jl995f ) to avo id spurious 



entrop y generation in regions with strong shear flows. The viscous term Hi,- was compu ted as in 
(I1992L Eqs. 4.1 and 4.2), and later corrected by a factor /y as defined in Balsara ( 19951 see also 
120051 ): 



Monaghan 



Monaghan 



TTCorr 

ij 



where 



fi 



^^Jh = H„0.5 (/, 



IV -t), 



(A2) 



(A3) 



V • w,| -f |V X -h0.0001c,//ii 

with V -Vi and V x the divergence and the curl of the velocity field, which we compute following Monaghanl 
(jl992l Eqs. 2.10 and 2.12), and Ci the sound velocity of particle i. 
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A. 3. Anisotropic interpolation kernel 

Although the use of artificial viscosity is not optimal in comparison with Riemann solvers, the use of 
anisotropic interpolation kernels greatly impr oves the abili ty of SPH cod es to deal with shock waves. We have 
adapted a 3D ellipsoidal kernel ([Shapiro et al. 1996: Ow en et al. I ll998h to our problem. When a detonation 



front is detected, the (otherwise spherically symmetric) interpolating kernel takes the shape of two half oblate 
spheroids joined at the equator (see Fig. [TO)) , sharing an axis of rotational symmetry that coincides with the 
direction of the local pressure gradient (computed with the standard spherically symmetric kernel) . The size 
of each semi-minor axis, 2ft,+ (2ft,_), is determined in three steps. A first value is obtained through: 

h±^Tn&-^[5i—^^,52-^r—], (A4) 



|vp±|' |vr±| 

where VP+ and Vr+ (VP_ and Vr_) are the pressure and temperature gradients at the position of particle 
i computed using only particles located in the forward (backward) direction of VP^ , and 5i and 82 are two 
coefficients of order unity that we take here as 5i — 82 — 0.13. The aim of the differentiation between 
both directions of VP^ is that the presence of a discontinuity in front (back) of a particle does not limit its 
capability to contribute to the determination of the smoothed quantities in the backward (forward) direction. 
In a second step we require that the flow be sonic by recomputing h± through: 

/i± = /i± max ^1,0.01— ^ , (A5) 

where c; and Vi are the sound velocity and the velocity of particle i, respectively. The last equation together 
with the pressure term in Eq. IA4I discard the activation of the anisotropic kernel in deflagration waves. 
Finally, the value of h± is constrained to be in the range h > h± > O.lh. In regions where there are no 
shocks the kernel reduces to the standard spherically symmetric kernel. In this way the spatial resolution at 
the detonation front improved, on average, by a factor '--^ 2 — 3. 

To compute the anisotropic kernel PK,^"' (u) at the position of a neighbor particle j (see Fig. [191) '^^ use the 



cubic spline kernel ([Monaghan fc Lattanzidll991l ) as a function of the ratio u between the particles distance 



rij and the size of the kernel in the direction of fij: u — 2rij/d, with: 



d^2 (A6) 
W /i± + (/i^ - /il) cos2 a 



and a the angle between fij and VP. The anisotropic kernel is finally normalized by a factor 

2n 



(A7) 



so that 



W^7dr=l. (A8) 



The ability of our device to follow the propagation of detonation waves was demonstrated in 
Garcia-Senz et al. ( 19991 Fig. 1), where more details of the numerical method are given, and in Figs. [31 
and [H in the present paper. A problem with ellipsoidal kernels is that they do not guarantee conservation of 
energy and angular momentum. In fact, during the detonation phase the violation of energy conservation 
reached ~ 5% of the final kinetic energy, while the net angular mom entum at the end of the simulation was 



negligibly small (for reference figures see Table 1 in lBravo Sz Garcfa-Senai2008i) 
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Table 1 

sph 3d delayed detonation models compared with other works 



Model 


Afdoti 
(Mo) 


K 

(10=1 ergs) 


M (^'^Ni) 
(Mq) 


M (NSE) 
(Mq) 


M (C + O) 

(Me) 


M (IME) 

(Mo) 


GCDl 


0.06 


1.51 


0.97 


1.21 


0.04 


0.10 


TURB7 


0.40 


1.49 


0.91 


1.19 


0.09 


0.07 


YlOO" 




1.52 


< 1.19 


1.19 


0.04 


0.13 


DD.005'' 


0.38 


1.52 


< 1.14 


1.14 


0.04 


0.22 



IPlewal | |2007|) 

I jRopke fc Niemeveii 1 120071) 
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Table 2 

Properties of the PRD models 



Model 


Jl-fdotl 


Mace" 




lldoto" 


K 






Amis" 




(Me) 


(Me) 


(1 000 km s~^) 


(1000 km s"^) 


{10" erg) 


(Mo) 


(erg s"^) 




PRD14 


0.14 


1.26 


14.0 


9.8 


1.22 


0.81 


2.4 X lO"*^ 


0.94 


PRD18'' 


0.18 


1.15 


11.0 


9.3 


1.05 


0.77 


2.0 X 10*^ 


1.32 


PRD26 


0.26 


1.05 


10.0 


8.0 


1.01 


0.64 


1.9 X lO"*^ 


0.93 


PRD29 


0.29 


0.83 















^Lagrangian mass of the accretion shock 

^Final ejecta velocity at the Lagrangian mass Mace 

'^Final ejecta velocity at the Lagrangian mass at which a confined detonation was ignited 
■^Maximum luminosity 

■^Decline in bolometric magnitude 15 days after bolometric maximum 

The nucleosynthesis of mode 1 PRD18, including the ^'^Ni yield, has been computed with a post-processing code 
a nd hence dif fers slightly from the figures given by the alpha- network built into the SPH code, which were provided 

in lBravo & Garcfa-Seri9 1I2006I) 
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Table 3 

Chemical composition of the e.tecta in solar masses 



Model 


PRD14 


PRD18 


PRD26 


C 


0.061 


0.101 


0.122 


o 


0.105 


0.149 


0.169 


Ne 


0.010 


0.015 


0.015 


Mg 


0.025 


0.027 


0.028 


Si 


0.118 


0.102 


0.138 


S 


0.055 


0.052 


0.068 


Ar 


0.011 


0.011 


0.014 


Ca 


0.010 


0.011 


0.013 


Cr 


0.005 


0.007 


0.006 


Mn 


0.001 


0.002 


0.001 


Fe 


0.878 


0.815 


0.720 


Ni 


0.054 


0.046 


0.051 
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Fig. 1. — Evolution of model PRD14 through 
the deflagration, pulsation, and detonation phases. 
Shown are the evolution of the total masses of fuel, 
intermediate-mass, and Fe-peak elements (left axis 
scale), and the central density (right axis scale) 




Fig. 2. — Evolution of temperature (color map, in 
10^ K) and fuel mass fraction (contours) in an oc- 
tant of model PRD14 during the deflagration and 
pulsation phases. The C-l-0 mass fraction con- 
tours shown are X = 0.3 (black), 0.5 (dark gray), 
and 0.8 (light gray). Axes units are thousands of 
km. Times shown are (from left to right and top 
to bottom): t = 0, 0.5, 1.0, and 2.0 s. Note that 
the size of the axes changes: 2 000 km in the top 
row images, and 6 000 km in the bottom row. The 
white dwarf radii at the times shown are: 1 900, 
2 300, 6 500, and 16 600 km 
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Fig. 3. — Detonation ignition in model PRD14. 
Temperature maps and fuel mass fraction con- 
tours are shown in a slice around the location of 
detonation formation: at t = 4.0 s, shortly af- 
ter ignition (top left image), and 0.1 s later (top 
right image). The bottom left image is a recon- 
struction of the hydrostatic force ( '^y- , color map 
in units of 10^*^ cm s^^) with density contours 
(log/9(g cm^"^) = 8, 7.5, 7, and 6.5 from the center 
outwards) at i = 4.1 s. The bottom right image 
shows in white the region in which a detonation 
wave has been detected by the detonation capture 
algorithm at t = 4.1 s. Axes units, temperature 
scale and fuel mass fractions represented are the 
same as in Fig. [21 The slice is a zoom on the yz 
plane of the octants shown in Figs. [2] and [5] 




Fig. 4. — Maximum resolution of the SPH cal- 
culation at i = 4.1 s in the same slice shown in 
Fig. m The color map represents the smoothing 
length, h, in km. The white points give the ac- 
tual position of the SPH particles for which the 
detonation capture algorithm has been activated 
and the anisotropic smoothing kernel is at work 
with its semiminor axis aligned with the pressure 
gradient 
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Fig. 5. — Dynamical evolution of model PRD14 
shortly before and during detonation propagation. 
Each thin line shows the radius of a spherical layer 
enclosing a constant mass. The thick lines show 
the inner and outer boundaries of the detonated 
layers, that are defined for the purposes of the 
present figure as those layers in which the mass 
of fuel processed by the detonation wave is larger 
than the mass of remaining fuel 



Fig. 6. — Evolution of temperature (color map, in 
units of 10^ K) and fuel mass fraction (contours) in 
an octant of model PRD14 during the detonation 
phase. The C+0 mass fraction contours shown 
are X — 0.3 (black), 0.5 (dark gray), and 0.8 (light 
gray). Axes units are thousands of km. Times 
shown are (from left to right and top to bottom): 
t = 3.8, 4.1, 4.4, and 4.7 s 
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Fig. 7. — Final angle-averaged density profile as a 
function of radius for model PRD14. The distri- 
bution of points resulting from the SPH calcula- 
tion has been divided into 1 000 radial shells. The 
mean density and the dispersion within each shell 
are shown 




1000 ''''' ' ' ' ' ' ' ' ' '' ' — ^ — ' — — 

5 10 15 20 25 

Velocity (1,000 km s-') 

Fig. 8. — Final angle-averaged density profiles of 
models PRD14 at i = 34 s (dotted hne), PRD18 
at 1.5 h (sohd hne), and PRD26 at 20 s (dashed 
line). The density is shown here factorized by 
in order to make the densities comparable 




Fig. 9. — Rendering of the final distribution of 
^^Ni in model PRD14. The isosurface defined by 
X(^^Ni) = 0.5 together with a volume rendering 
of the ejecta density are shown 




Fig. 10. — Final chemical composition in an oc- 
tant of model PRD14. It is interesting to note that 
the chemical structure at the end of the pulsating 
phase (Fig. 3 in Paper I) left no imprint on the 
final distribution of elements. The detonation led 
to a quite symmetric distribution with some de- 
gree of chemical stratification, with C-0 confined 
to the outskirts of the ejecta, Si-S dominant at 
intermediate distances, and Fe-group nuclei con- 
centrated in the inner ejecta but also present in 
clumps distributed through the outer ejecta. The 
axis units are thousands of km 
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Velocity (1 000 km s-') 



Fig. 11. — Final angle- averaged chemical profiles 
of model PRD14 as a function of the velocity. The 
distribution of points resulting from the SPH cal- 
culation has been divided into 1 000 radial shells. 
The mean mass fraction and the dispersion within 
each shell are shown. The fuel (C-l-0) mass frac- 
tion is painted in blue, that of intermediate-mass 
elements in green, and that of Fe-group elements 
in red 
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Velocity (1 000 km s"') 



Fig. 12.— Same as Fig.[TIl but for model PRD18 




5 10 15 20 25 

Velocity (1 000 km s"') 



Fig. 13.~ Same as Fig.ITTl but for model PRD26 



T 




Velocity (1 ,000 km/s) 

Fig. 14. — Abundances at t = oo in the ejecta of 
model PRD14 as a function of the velocity. The 
abundances are angular averaged in 130 spherical 
shells 
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Lagrangian mass coordinate (M^) , 

Am, 5 (bo. 



Fig. 15. — Chemical profile of the delayed det- 
onation model TURB7 just after the end of the 
detonation phase (points) and in the homologous 
expansion phase (lines: t = 62 s). Carbon and 
oxygen are represented by blue triangles and the 
dotted line, intermediate-mass elements by green 
squares and the dashed line, and Fe-peak elements 
by red circles and the solid line 



Fig. 16.— Mass of ^^Ni synthesized in SNIa 
explosions vs decrease in bolometric magnitude 
15 days after maximum. Starred symbols are 
data for a sample o f well observed SNIa from 
Stritzinger et al.l (|20061 full squares belong to the 
Pulsating Reverse Detonation models, while full 
triangles represent models GCDl and TURB7 
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Projected velocity (1,000 km s ^ 



Fig. 17. — Column density of Fe-group elements 
for model PRD14 at day 25 after thermal runaway. 
The column density (in g-cm~^) of Fe-group ele- 
ments between the photosphere and the observer 
is color-coded for an arbitrary line of sight. White 
and red colors mark places for which the mate- 
rial above the photosphere is nearly pure Fe, while 
blue and green indicate the contrary: that Fe is 
scarce. 



Fig. 18. — Distribution of the mass of carbon be- 
tween the photosphere and the observer as a func- 
tion of the velocity in the line of sight for model 
PRD14 at day 13 after thermal runaway. The dif- 
ferent curves belong to 5 arbitrary lines of sight. 
For this calculation the ejecta has been mapped to 
a three-dimensional grid in spherical coordinates, 
and the location of the photosphere has been esti- 
mated as for Fig. [T7] (see text) 



26 



Fig. 19. — Geometry of the ellipsoidal kernel of 
particle i. The semi-major axis along a plane 
normal to the pressure gradient is the smoothing 
length h, while the semi-minor axis is aligned with 
Vp. In our implementation we use two different 
values of the semi-minor axis, one in each direction 
of the pressure gradient: and /i_ (see text for 
details). The kernel at the position of a neighbor 
particle j, depends on the angle, o, between the 
relative position vector r and the direction of the 
semi-minor axis, and on the ratio u of the parti- 
cles distance to the size of the ellipsoid, d, in the 
direction of r: u = 2r/d 



